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>"^ Aims. The Spectral and Photometric Imaging REceiver (SPIRE) onboard the Herschel space telescope has provided confusion limited maps of 
\^J deep fields at 250/im, 350 /an, and 500 |j.m, as part of the Herschel Multi-tiered Extragalactic Survey (HerMES). Unfortunately, due to confusion, 
only a small fraction of the cosmic infrared background (CIB) can be resolved into individually-detected sources. Our goal is to produce deep 
galaxy number counts and redshift distributions below the confusion limit at SPIRE wavelengths (~20 mjy), which we then use to place strong 
constraints on the origins of the cosmic infrared background and on models of galaxy evolution. 
Q Methods. We individually extracted the bright SPIRE sources (>20mJy) in the COSMOS field with a method using the positions, the flux 
I— i densities, and the redshifts of the 24 jim sources as a prior, and derived the number counts and redshift distributions of the bright SPIRE sources. 
TJ* For fainter SPIRE sources (<20 mjy), we reconstructed the number counts and the redshift distribution below the confusion limit using the deep 
C^J 24 |j.m catalogs associated with photometric redshift and information provided by the stacking of these sources into the deep SPIRE maps of the 
1 ' GOODS-N and COSMOS fields. Finally, by integrating all these counts, we studied the contribution of the galaxies to the CIB as a function of 

their flux density and redshift. 

~~^ Results. Through stacking, we managed to reconstruct the source counts per redshift slice down to ~2 mjy in the three SPIRE bands, which lies 

^ about a factor 10 below the 5tr confusion limit. Our measurements place tight constraints on source population models. None of the pre-existing 

1*1 models are able to reproduce our results at better than 3-cr. Finally, we extrapolate our counts to zero flux density in order to derive an estimate 

CN of the total contribution of galaxies to the CIB, finding lO.l^'nWm -2 sr" 1 , 6.5^|-gnWnT 2 sr _1 , and 2.8^-gnWnT 2 sr~' at 250 |im, 350 \xm, and 

^^ 500 urn, respectively. These values agree well with FIRAS absolute measurements, suggesting our number counts and their extrapolation are 

sufficient to explain the CIB. We find that half of the CIB is emitted at z =1.04, 1.20, and 1.25, respectively. Finally, combining our results with 

#yv other works, we estimate the energy budget contained in the CIB between 8 p.m and 1000 (im: 26^ nW nT 2 sr~' . 

Key words. Cosmology: observations - Cosmology: diffuse radiation - Galaxies: statistics - Galaxies: photometry - Submillimeter: galaxies - 
Submillimeter: diffuse background 
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1. Introduction (~15%) are due to accretion processes (Alexander et al. 2005; 

Jauzac et al. 2011). The CIB thus primarily gives a budget of 

C3 About half of the relic energy arising from the emission of infrared photons emitted by star-formation processes, 
galaxies, which we refer to as the Extragalactic Background 

Light (EBL), is contained in the cosmic infrared background More recently, deep number counts (flux density dis- 

(CIB), which lies between 8 urn and 1000 urn, and peaks at tributions of infrared sources) have been measured in the 

around 150 urn (Hauser & Dwek 2001; Dole et al. 2006). The mid- and far-infrared domain, thanks to the sensitivity of the 

first absolute measurements of the CIB were performed in Spitzer (Werner et al. 2004) and Herschel 1 (Pilbratt et al. 

the nineties with the Far-Infrared Absolute Spectrophotometer 2010) space telescopes. They exhibit power-law behavior at 

(FIRAS; Puget et al. 1996; Fixsen et al. 1998; Lagache the faint end, which can be extrapolated to zero flux density 

et al. 1999) and the Diffuse Infrared Background Experiment in order to estimate the contribution of all the galaxies to the 

(DIRBE; Hauser et al. 1998) onboard the COsmic Background 



Explorer (COBE). The far-infrared emission from galaxies is i Herschel is an ES A space observatory with science instruments pro- 
mainly due to dust heated by ultraviolet photons re-radiate in vided by European-led Principle Investigator consortia and with impor- 
the infrared. A small fraction of these far-infrared emission tant participation from NASA. 
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CIB (e.g. Papovich et al. 2004 at 24 pm with Spitzer/MlPS, 
Bethermin et al. 2010a at 24 pm, 70 pm, and 160 pm with 
Spitze r/MIPS, Berta et al. 201 1 at 70 pm, 100 pm, and 160 pm 
with Herschel/PACS). These estimations of the total CIB 
agree with the absolute measurements performed by COBE, 
suggesting the CIB is now explained shortward of 160 pm. At 
longer wavelengths, due to source confusion (Dole et al. 2003; 
Nguyen et al. 2010), the Herschel/SPIRE instrument (Griffin 
et al. 2010) can directly resolve only 20%, 12%, and 6% of the 
CIB at 250 \xm, 350 |_im, and 500 [am, respectively (Oliver et al. 
2010b). Due to the limited depth of these confusion-limited 
observations, the break and the power-law behavior of the 
counts at faint flux density cannot be seen. It is thus necessary 
to use statistical tools like P(D) analysis 2 (Condon 1974; 
Patanchon et al. 2009) or stacking (Dole et al. 2006; Marsden 
et al. 2009) to study the origins of the sub-mm part of the CIB. 

Using a P(D) analysis, Patanchon et al. (2009) pro- 
duced deep counts from the Balloon-borne Large- Aperture 
Submillimeter Telescope (BLAST, Pascale et al. 2008; Devlin 
et al. 2009) data. They were only able to constrain one data 
point below 100 mJy at 250 [J.m, and were not sensitive to more 
subtle features in the shape of the counts. Using a stacking 
analysis, Bethermin et al. (2010b) managed to detect the peak of 
the Euclidian-normalized counts at 250 \xm in the BLAST data, 
but not at longer wavelengths. Using a P(D) analysis on SPIRE 
data, Glenn et al. (2010) managed to clearly detect this peak 
at 250 \xm and 350 \im, but not at 500 [am. In all these cases, 
the uncertainties are too large to reliably detect a power-law 
behavior at the faint end. 

A stacking analysis of SPIRE data similar to that performed 
on BLAST data by Bethermin et al. (2010b) could significantly 
reduce the uncertainties and provide more precise information 
on the sources which make up the CIB. Le Floc'h et al. (2009) 
and Berta et al. (2011) also showed that counts per redshift 
slice are strong constraints for galaxy evolution models (e.g. Le 
Borgne et al. 2009, Valiante et al. 2009, Marsden et al. 2011, 
Bethermin et al. 2011, Gruppioni et al. 2011, Rahmati & van 
der Werf 2011). Lastly, unlike a P(D) analysis, stacking allows 
us to measure directly the counts in redshift slices, but requires 
a prior catalog. Thus, here we perform a stacking analysis in the 
SPIRE bands, in the COSMOS and GOODS -N fields to produce 
deep counts per redshift slice in SPIRE bands, combining the 
Herschel Multi-tiered Extragalactic Survey (HerMES) 3 data 
(Oliver et al. 201 1) and the ancillary data. 

The paper is organized as follows. In Sect. 2, we present the 
different data sets used in our analysis. We then introduce the 
method used to measure the number counts of resolved sources 
(Sect. 3) and another method based on stacking to reconstruct the 
number counts below the confusion limit (Sect. 4). In Sect. 5, 
we detail the estimation of the statistical uncertainties. Sect. 6 
presents a end-to-end simulation used to check the accuracy of 
our method. In Sect 7, we interpret our number counts and com- 
pare them with previous measurements and models of galaxy 
evolution. The same thing is done in Sect. 8 for the redshift dis- 
tributions. In Sect. 9, we derive constraints on the CIB level and 
its redshift distribution from our number counts. We finally dis- 
cuss our results (Sect. 10) and conclude (Sect. 11). 



2 P(D) analysis is a statistical method used to estimate the number 
counts in a field from the pixel histogram of an extragalactic map. 

3 hermes.sussex.ac.uk 



2. Data 

2.1. SPIRE maps at 250 \im, 350 \im and 500 p.m 

The SPIRE instrument (Griffin et al. 2010) onboard the Herschel 
Space Observatory (Pilbratt et al. 2010) observed the COSMOS 
field as part of the Herschel Multi-tiered Extragalactic Survey 
(HerMES) program Oliver et al. (2011). The maps were built 
using an iterative map-making technique (Levenson et al. 2010). 
The full width at half maximum (FWHM) of the SPIRE beam 
(Swinyard et al. 2010) is 18.1", 24.9", and 36.6" at 250 \im, 
350 \im, and 500 |am, respectively. The typical instrumental 
noise is 1.6, 1.3, and 1.9mJybeam^ in COSMOS (1.6, 1.3, and 
2.0mJybeam _1 in GOODS-N) and the 1-cr confusion noise is 
5.8, 6.3 and 6.8mJybeairT 1 in the three wavebands (Nguyen 
et al. 2010). The maps are thus confusion limited. The absolute 
calibration uncertainties in point sources are estimated to be 
7% (Swinyard et al. 2010, updated in the SPIRE Observers' 
Manual 4 ). 



2.2. Ancillary data In COSMOS 

Deep 24 ixm imaging of the COSMOS field was performed 
by the Spitzer Space Telescope (S-COSMOS, Sanders et al. 
2007). The associated catalog reaches 90% completeness at 
80 iiJy (Le Floc'h et al. 2009). This catalog was matched with 
photometric redshifts of Ilbert et al. (2009). Due to the high 
density of optical sources compared with the size of the MIPS 
beam, the cross-identification can be ambiguous in many cases. 
An intermediate matching was thus performed with the K and 
IRAC bands where the source density is smaller, which helps 
to discriminate between several optical counterparts in a MIPS 
beam (Le Floc'h et al. 2009). 

In this paper, we use an updated version of the photometric 
redshift catalog of Ilbert et al. (2009) (vl.8). This version uses 
new deep //-band data. However, this catalog is not optimized 
for AGN. For the sources detected by XMM-Newton we instead 
use the photometric redshifts of Salvato et al. (2009), estimated 
with a technique specific to AGN. In addition, 10 000 sources 
have spectroscopic redshifts provided by the S-COSMOS team 
Lilly et al. (2007), which, where available, are used instead 
of the photometric redshifts. Details of the updated COSMOS 
5 24 + z catalog will be given in Le Floc'h et al. (in prep.). In 
this new version, 96% of the 27 811 S 24 > 80/xly sources have 
redshifts (9.7% of them are spectroscopic). 



2.3. Ancillary data in GOODS-N 

In the GOODS-N field, we use the 24 11m catalog of Magnelli 
et al. (2011). This catalog was built using the IRAC catalog at 
3.6 [am as a prior, and has an estimated 3-cr depth of 20 piy, 
but at this depth the completeness is only ~50%. The stacking 
of an incomplete catalog can bias the results (Bethermin et al. 
2010b; Heinis et al. in prep., Vieira et al. in prep.). According 
to simulations, cutting at 80% completeness results in smaller 
bias. We thus cut the catalog at 30/iJy (80% completeness limit) 
to have a more complete and reliable sample. These sources 
were matched with the photometric redshifts of Eales et al. 
(2010) (97.4% of the 2791 24pm sources are associated with a 
redshift). The data fusion of these two catalogs will be explained 



http://herschel.esac.esa.int/Docs/SPIRE/pdf/spire_om.pdf 
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in more detail in Vaccari et al. (in prep.). 



3. Measuring the statistical properties of the 
resolved sources 

In order to build counts per redshift slice and redshift distribu- 
tions of the sources selected by their SPIRE flux densities, we 
require catalogs containing SPIRE flux densities and redshifts. 
The redshift catalogs are built from optical and near-infrared 
catalogs. We start from catalogs of the 24 [am sources which 
have optical counterparts and thus photometric redshifts. Due to 
the large beam of SPIRE, it is not trivial to identify the MIPS 
24 (j.m counterpart for a given SPIRE source. To avoid this 
problem, we directly measure the SPIRE flux denisty of the 
24 [im sources in the maps by PSF-fitting assuming a known 
position (Bethermin et al. 2010b; Chapin et al. 201 1; Roseboom 
et al. 2010). Sect. 3.3 and 7 discuss the relevance of this choice 
of prior. 

The GOODS-N field, being much smaller than COSMOS, 
has little impact on the statistical uncertainties (using GOODS- 
N+COSMOS reduces the uncertainties of 0.2% compared 
to COSMOS only). The inclusion of GOODS-N introduces 
heterogeneity to the 24 um catalogs, which are built using 
different methods between the two fields. For this reason, we 
used only the COSMOS field in the following section. 



3.1. Source extraction 

We use the fastphot PSF-fitting routine, described in Bethermin 
et al. (2010b), to fit the following model to the SPIRE map: 
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where m is the map, N s the number of sources, S^ the SPIRE 
flux density of the £-th source, b Xtt y k a point spread function 
(PSF) centered on the position of the k-th source {xk,yk), and /j. 
a constant background. The catalog of positions used as input to 
fastphot is discussed below. The free parameters fit by fastphot 
are the SPIRE fluxes of sources in the prior list Sk and the 
level of the constant background fi. We used the PSF based on 
the Neptune scan from Glenn et al. (2010) 5 . The map is not 
fit in one pass, but split into 100x100 pixel regions (the pixel 
sizes are 6.0", 8.3", and 12" at 250 (i.m, 350 |im, and 500 [J.m, 
respectively). Each region was fit independently. To limit the 
edge effects, we also fit simultaneously an additional region 
of 20 pixels around each 100x100 region. The positions of the 
sources in both the central and additional regions are used by 
fastphot, but we keep in the final catalog only the photometry 
of the sources in the central region. The signal at 20 pixels (~6 
times the FWHM) from the center of a source is negligible. A 
source outside of the additional region cannot thus significantly 
affect the photometry in the central region. 

The fastphot routine suffers some instabilities when two 
sources are too close to one another. We thus do not use the 
position of all the 24 yum sources in fastphot. For several 



5 Beam data are also available from the Herschel Science Centre at 
ftp://ftp.sciops.esa.int/pub/hsc-calibration/SPIRE/PHOT/Beams 



redshift and 24 [am flux density slices, we estimated the mean 
color by stacking (see Sect. 4.1). We then use these mean 
colors to estimate the flux density of each source in the 
SPIRE bands. A 24 inn source is included in the position list 
of fastphot only if it has the highest estimated SPIRE flux 
density in a 0.5xFWHM radius. This process was therefore 
performed independently in each band. Some sources with 
unusually high sub-mm/mid-infrared colors could be missed by 
this method, but there are few objects of this kind (see Sect. 3.3). 

To avoid unphysical negative flux densities for faint sources 
lying on negative fluctuations of the noise, we run fastphot 
iteratively, removing from the position list the sources with 
negative flux densities at each iteration. Removing a source 
from the input catalog is equivalent to assuming its flux density 
is zero, which is the most probable value in this case. 



3.2. Estimating photometric noise 

We estimate the photometric noise from the standard deviation 
of the fastphot residual map, finding the values 4.6, 5.5 and 
5.1 mJy at 250 inn, 350 jam, and 500 [J.m, respectively. These 
values are about 20% lower than the combination the l-<x 
confusion noise measured by Nguyen et al. (2010) and the 
instrumental noise (6.0, 6.4 and 7.0 mJy). Our method is thus 
more efficient than a naive blind extraction. We chose to cut 
our statistical analysis at 20 mJy in the three SPIRE passbands, 
which corresponds to about A-cr. 

In order to cross-check our estimate of the photometric 
noise, we inject 200 artificial point sources in the real SPIRE 
map and add them in the input position list of fastphot. We add 
a random shift drawn from a 2D Gaussian with <x = 2" to the 
source position in order to simulate the astrometric uncertainties 
of the real catalog. We then rerun fastphot and compare the 
input and output flux densities. Fig. 1 shows the histogram of 
the difference between the recovered and input flux densities. 
We found a l-cr photometric noise of 3.9, 5.2 and 5.1 mJy at 
250 (o.m, 350 [J.m, and 500 \xm, respectively. The values are 
similar to those estimated from the residual map. Comparing 
the two sets of values, we can estimate an error of about 20% on 
the photometric noise. 



3.3. 24 |xm Dropouts 

A fundamental limitation of our model is that it is not sensitive 
to any population of sources that are faint at 24 microns but 
bright in the SPIRE passbands (24 fj.m dropouts). No such 
population is known or theoretically predicted, except possibly 
at very high redshifts, but the possibility remains that nature has 
been more inventive than we have. In this section we attempt to 
test whether there is any evidence for such sources, and do not 
find any. 

First, we study the residual SPIRE maps after removing 
all of the sources extracted with FASTPHOT to estimate the 
number of sources missed by the extraction. The density of 
remaining sources is quite small so that confusion is not a 
problem. We then search for additional sources by looking for 
peaks in the beam smoothed residual map to a depth of 20 mJy, 



Bethermin et al.: Deep number counts at 250/mi, 350/mi and 500 u.m and CIB build-up. 



40 


; ' ' ' ; 




250 fim 


30 


- 350/mi 






- 




500 /mi 


l 






— I 


* 20 














-_ 


10 












~b 




- 










; 




:...F=W. 








=l , 




; 























,..l, 



-20 



-10 



10 



20 



S -S [mJy] 

out in 



Figure 1. Simulation of the photometric uncertainties: histogram of the 
difference between the recovered and the input flux density of the ar- 
tificial sources injected in the real SPIRE map and re-extracted with 

FASTPHOT. 



and find that the fraction 6 of such possible sources is only 1.3%, 
0.7%, and 0.6% at 250 /mi, 350 /mi, and 500 /mi, respectively. 
Next, we have compared our prior catalog with the source list 
from the blindly-extracted HerMES catalog (Smith et al. 2012), 
which is limited to sources brighter than 20 mJy. The fraction 
of sources in the blind catalog without counterpart in the prior 
catalog strongly depends on the choice of matching radius. For 
a narrow radius of 0.5 FWHM, we obtained 3.2%, 2.6%, and 
0.6% whereas for a large radius of 1 FWHM, we obtained 0.6%, 
0.3%, and 0.0%, at 250 /mi, 350 /mi, and 500 /urn, respectively. 
However, with the narrow radius, we miss some sources due to 
astrometric uncertainties, and with the large radius, we possible 
have a contamination by neighboring sources. It is expected 
that the fraction of dropouts decreases with the flux density. 
Nevertheless, this behavior is hard to constrain because of the 
small number of bright sources. Note however that the fraction 
of dropouts at the flux density limit is very close to the values 
obtained for the full sample because of the steep slope of the 
counts. Thus, we conclude that our catalog, based on the 24 |mi 
prior, is very close to complete above 20 mJy. From a blind 
extracted catalog of H-GOODS data, Magdis et al. (2011) 
estimated the fraction of 24 //m dropouts for H-GOODS fields 
and shallower fields. They predict a dropout faction smaller than 
2% in the COSMOS field. 

Finally, we can compare our number counts measurement 
to other analyses that did not make use of a 24 /mi prior. 
We find good agreement with the blind extractions of Oliver 
et al. (2010b) and Clements et al. (2010), but note that these 
analyses required significant model corrections for Eddington 
bias and confusion. We also find good agreement with the P(D) 
analysis of the SPIRE maps by Glenn et al. (2010), which, 
by construction, is not affected by either issue. We take these 
comparisons as a strong indication that we have not missed a 
statistically significant population, at least in terms of the red- 



The fraction of dropouts is defined as 



N t l+N p 



. where N d is the num- 



ber of sources brighter than 20 mJy found in the residual map and N p 
the number of sources brighter than 20 mJy extracted by fastphot in the 
normal map. 
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Figure 2. Effect of the photometric uncertainties in the flux density 
distribution at 500 \im. Black dashed line: distribution of the flux den- 
sity measured at the position of the 24 p.m sources. Red solid line: the 
same distribution after adding a 5. 1 mJy random Gaussian noise to each 
measured flux density. Due to photometric noise, some sources have a 
negative flux density (put to zero in our iterative algorithm) and are not 
represented here. Black dotted line: flux density cut used in our analysis 
(20 mJy). 



shift integrated number counts. However, we must acknowledge 
that the fraction of dropouts could evolve with redshift, and in 
particular the high redshift bins may be less complete. 



3.4. Correction of the biases 

Intuitively, the simplest way to compute the source counts is 
to measure the number of sources in a flux density bin and 
divide it by the width of the bin and the surface area of the field. 
However, due to photometric noise, this estimate is biased. In 
fact, for a prior-based extraction, we do not have a flux boosting 
phenomena (which appears at low signal to noise ratio for 
blind extraction, because the completeness is higher for sources 
lying on positive fluctuations of the noise, as discussed e.g. in 
Bethermin et al. 2010b), but another statistical effect, Eddington 
bias, biases the counts measurement, as illustrated by Fig. 2. 
The black dashed line shows the distribution of the 500 /mi 
flux densities measured at the prior positions. We will assume 
this distribution is close to the real one, and will somewhat 
arbitrarily refer to it as initial distribution. The red line shows 
the same distribution, but adding a 1 -<x 5.1 mJy Gaussian error 
on the flux density of each source, called measured distribution. 
At bright flux density (S500 > 20mJy), we can observe an 
excess in the measured distribution compared to the initial one. 

To correct this bias, we use a Monte Carlo (MC) method as 
in Bethermin et al. (2010b). We compute 1000 realizations of the 
bias in each flux density (regular in logarithm from 20 mJy) and 
redshift (0 < z < 0.5, 0.5 < z < 1, 1 < z < 2, and z > 2) bin, and 
use them to compute the mean correction and its uncertainty: 

- We start from the measured distribution of the 250 /mi, 
350 /mi, or 500 //m flux density of the 24 //m sources in the 
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prior list of fastphot in a given redshift bin, and assume it 
is close to the initial distribution. This last hypothesis is a 
significant approximation, but the selection function of the 
procedure used to construct the prior catalog is too complex 
to be modeled without introducing strong assumptions about 
galaxy evolution. 

- We draw with replacement N sources in the initial sample, 
where N is the number of sources in the initial sample. This 
bootstrap step is used to take into account the sample vari- 
ance on the initial flux density distribution. 

- We add a Gaussian random photometric noise to the flux 
density of each source. We use the values of the noise found 
in Sect. 3.2 plus a 20% systematic shift (different at each it- 
eration of the MC procedure), which takes into account the 
systematic uncertainty on the determination of the noise. 

- We compute the bias on the counts dividing the counts from 
the drawn sample before and after adding the photometric 
noise. 

Table B.l shows the corrective factor in various flux density 
and redshift bins and at various wavelengths. This correction 
can reach 40% in the fainter flux density bins and decreases at 
brighter flux densities. 
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Figure 3. Mean flux density measured by stacking as a function of 
wavelength. The various redshift bins are represented using various col- 
ors. Each panel corresponds to each 24 /mi flux density bins. The error 
bars are estimated with a bootstrap method. 



Similar corrections are applied to the redshift distributions. 
However, in addition, we apply a random error to the redshift 
based on the uncertainties provided in the photo-z catalog (but 
without taking into account the catastrophic outliers) during the 
MC procedure. In this case, there is only one flux density bin 
(>20 mJy). 



4. Measuring the statistical properties of sources 
below the confusion limit 

Due to source confusion, SPIRE cannot resolve the bulk of 
the CIB into individual sources (Nguyen et al. 2010; Oliver 
et al. 2010b). Nevertheless, about 80% of the CIB is resolved 
at 24 u.m (Papovich et al. 2004; Bethermin et al. 2010a). We 
thus perform a stacking analysis using the 24 urn prior to probe 
fainter populations and resolve a larger fraction of the sub-mm 
CIB. 



4.1. Stacking method 

Stacking is a statistical method which allows us to measure the 
mean flux density of a population of sources selected at another 
wavelength, but which are too faint to be detected individually 
at the working wavelength. Several methods can be used (e.g. 
Dole et al. 2006, Marsden et al. 2009, see the discussion in 
Vieira et al. in prep.). We use the following method (also used 
in Vieira et al. in prep.): we first subtract the mean of the SPIRE 
map in the region covered by the 24 urn observations. We then 
compute the mean signal in pixels which has a source centered 
on them. This provides the mean flux density of the population, 
because the SPIRE maps are in Jy beairT 1 . Vieira et al. (in prep.) 
showed that this method is more accurate in a confusion-limited 
case than PSF-fitting on a stacked image. The uncertainties are 
estimated using a bootstrap method (Jauzac et al. 201 1). 

Due to the large number of sources in COSMOS, we can 
split our 24 urn sample into eight redshift bins (0 < z < 0.25, 
0.25 < z < 0.5, 0.5 < z < 0.75, 0.75 < z < 1, 1 < z < 



1.5, 1.5 < z < 2, 2 < z < 3, and z > 3) and logarithmic flux 
density slices (80/iJy< 5 24 < 172/Jy, 172/dy< S24 < 371 /jJy, 
371 /dy< S 2 4 < 800 /Jy, and 800 //Jy< S 2 a < 1723 //Jy). In 
GOODS-N, we use the same redshift slices, but a single flux 
density slice (30//Jy< S24 < 80/Jy). This choice of the number 
of bins was done to have a compromise between a fine grids in 
24/im flux density and redshift, but also a reasonable number 
of sources to stack in each bins to obtain a good signal-to-noise 
ratio. We stack the sources in each bin to compute their mean 
flux density in the three SPIRE bands. Fig. 3 shows the mean 
flux density as a function of wavelength, which, as expected, 
decreases rapidly in low redshift bins and peaks between 350 /mi 
and 500 /mi in the z > 3 bin. The mean color in each bin is 
computed by dividing the mean SPIRE flux density by the mean 
24 //m flux density. 

4.2. Scatter of the photometric properties of the stacked 
populations 

The uncertainties given by the bootstrap method, cn, oot , are 



°"boot 



4 



cr 2 , + cr 2 f + cr 2 

instr cont P°P 



Va^ 



(2) 



where cr, ns/r is the instrumental noise, cr con f the confusion noise, 
cTp p is the standard deviation of the flux density of the popu- 
lation, and Attack the number of stacked sources. The quantity 



4 



o"r , + cr e can be estimated from the standard deviation of 

instr conr 

the map, allowing us to deduce <x pop from our bootstrap analysis. 

While this formula is true for a Gaussian distribution, we 
note that the distribution of colors of the sources are probably 
better described by a log-normal distribution. Fig. 4 shows 
the distribution of the logarithm of the S 250/524 colors of the 
resolved sources (S 250 > 20 mJy) in the 1 < z < 1.5 redshift bin 
(this redshift bin was chosen because it has the larger number 
of sources). The confusion noise is also non Gaussian (Glenn 
et al. 2010). Nevertheless, due to central limit theorem, these 
distributions of the mean flux density tend to be Gaussian if 
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Figure 4. Black histogram: distribution of the logarithm of the 5 250 /S 24 
colors of the resolved sources (5 250 > 20mJy) in the l<z<1.5 redshift 
bin. Red line: fit of the histogram by a Gaussian. 
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Figure 5. Red histogram: Pixel histogram of the 250 jim SPIRE map in 
COSMOS. Blue histogram: Histogram of the mean signal in 100 pixels 
taken randomly in 100 000 realizations. Red and blue lines: Gaussian fit 
of the previous histograms. 



agrees with the value of 62% found for the resolved sources (see 
Fig. 4). 



4.3. Does the clustering of sources introduce a bias? 

The simplest stacking method assumes implicitly that the 
sources in the map are not clustered, but this has been shown 
to be unrealistic and must be accounted for (Bethermin et al. 
2010b; Viero et al. 2011; Penner et al. 2011). We have per- 
formed several tests in the COSMOS field to estimate the bias 
due to clustering, which we now describe. 



4.3.1 . Method A: convolution of the 24 urn map with the 
SPIRE beam 

A simple way to estimate the bias due to clustering is to 
convolve the 24 \xm map with a Gaussian kernel to obtain a 
24 i4.m map with a Gaussian PSF of the same FWHM as the 
SPIRE map Oliver et al. (2010a). To match resolutions, we use a 

We 



2 2 

""SPIRE _ °"mIPS- 



Gaussian kernel with beamsize enamel = J°~s 
measure the mean flux density of the 24 \im sources by stacking 
the 24 /im catalog on this convolved 24 (a.m map. The bias due 
to clustering is estimated by comparing the mean flux density 
measured by stacking with the mean flux density estimated from 
the 24 /mi catalog. We find biases of 5±2%, 11+2%, and 20+5% 
at 250 |im, 350 |J.m, and 500 \im, respectively. This method 
is equivalent to building a simulated map assuming a single 
color C for all the 24 /zm objects (including the ones below 
the detection limit at 24 yum), measuring the mean flux density 
of the selected population by stacking on this convolved map, 
and comparing it with the mean flux density coming from the 
catalog (< 5 24 > xC). The same color factor C is present in the 
mean stacked flux measured by stacking in the convolved map 
and in the mean flux coming from the catalog. It thus disappears 
when we compute the relative difference between these two 
quantities. We thus take C — 1 for simplicity. As expected, the 
bias due to clustering increases with the size of the beam. This 
estimate is exact only if the S spire/S24 color is constant, or if 
the properties of the angular clustering do not evolve with the 
color of the sources (and thus the redshift). These assumptions 
are not going to be exactly met, so that next we use another 
method to cross-check this estimate. 



a sufficient number of sources are stacked. Fig. 5 illustrates 
this property. The red histogram is the pixel histogram of the 
250 //m SPIRE map. It is not Gaussian, because the confusion 
noise is not. The blue one is the distribution of the mean signal 
in 100 pixels taken randomly in 100000 realizations (typically 
the effect of the instrumental and confusion noise on a stack of 
100 sources). This histogram is much closer to a Gaussian. The 
same thing happens for the color scatter term. The Gaussian 
approximation is thus very relevant here. 

The scatter on the S spire IS 24 color can be estimated by 
dividing cr pop by the mean flux density of the population. 
We do not detect a significant evolution of this scatter with 
redshift, wavelength or 24 (j.m flux density. We use the median 
and the standard deviation of the values found in the differ- 
ent redshift and 24 |im flux density bins, and find a scatter 
o"coior = 0" P op/ < S spire > of 68±35% (in linear units). This 



4.3.2. Method B: simulation based on mean colors 
measured by stacking 

In Sect. 4.1, the S spire /S 24 mean color as a function of the 
24 i4.m flux density and redshift were measured by stacking. We 
use these mean colors and the scatter measured in Sect. 4.2 to 
generate mock SPIRE flux densities for the sources in the S 24 +z 
catalog, and then build a simulated map of the COSMOS field 
using the position given in the 5 24 + z catalog, the estimated 
SPIRE flux density, and the SPIRE PSF. Random Gaussian 
noise was added following the noise map of the real data. We 
then stacked all the 24 \±m sources, and compared the mean 
flux density measured by stacking in the simulated map and 
the mean flux density in the mock catalog. We find a bias of 
7.0±0.9%, 10.4± 0.7%, and 20.6+1.2% at 250 |xm, 350 ^m, 
and 500 \im, respectively, in agreement with the values pro- 
vided by method A. The main drawback of this method is that 
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Figure 6. Radial profile of the stacked image at 250 \im of all the 
27 811 524 > 80 /uJy sources in COSMOS. Black squares: measure- 
ments. Red solid line: best fit. Green dot-dashed line: contribution of 
the PSF. Blue dashed line: contribution of the clustering. The error bars 
are too small to be represented. The pixel size is 6". 

any bias due to sources undetected at 24 microns is not modeled. 

We have also stacked sub-samples selected in redshift and/or 
in 24 u.m flux density. The bias tends to slightly decrease with 
the redshift and the 24 u.m flux density cut. Nevertheless, this 
evolution is small (below 3%), and the significance is smaller 
than 3-cr, We thus chose to neglect it, and assume a single value 
for the bias due to the clustering. 



4.3.3. Method C: fitting the profile of the stacked image 

For method C, we follow the Dole et al. (2006) method to pro- 
duce our stacked images. In the presence of clustering, this im- 
age can be fit by the following function (Bethermin et al. 2010b; 
Henis et al. in prep.): 



M - axb +/?x 



b * w 



max(b * w) 



(3) 



where M is the stacked image, w the auto-correlation function 
(ACF), * the convolution product, and b the beam function. The 
PSF is normalized to unity at the center to match the per-beam 
normalization of the maps, a and (3 are free parameters in the fit. 
The results of the fit are plotted in Fig. 6. In order to estimate the 
uncertainties, the fit was performed on 1000 bootstrap samples. 
If we measure the photometry in the central pixel of the PSF, the 
bias due to clustering is/3/a. We found 7.7+0.5%, 10.3+0.8%, 
and 19.1 + 1.8% at 250 u.m, 350 urn, and 500 u.m, respectively. 
The uncertainty here is the standard deviation of the values found 
for the different bootstrap samples. As expected, we also found 
that a and /3 are significantly anti-correlated (correlation coef- 
ficients of -0.46, -0.56, and -0.57 at 250 u.m, 350 urn, and 
500 urn, respectively). As was the case for method B, we do not 
detect any significant evolution of the bias with redshift. 

4.3.4. Correction of the bias due to clustering 

Table 1 summarizes our estimates of the bias due to clustering. 
Our three methods give similar results. To correct for the effects 



250 
350 
500 



5±2% 7.0±0.9% 7.7±0.5% 
11±2% 10.4+0.7% 10.3+0.8% 
20+5% 20.6+1.2% 19.1+1.8% 



Table 1. Bias due to clustering as a function of the wavelength. These 
values are estimated with the methods presented in Sect. 4.3. 



of the clustering, we divide our measured mean flux densities 
by the mean values of 1.07, 1.10, and 1.20 at 250 u.m, 350 u.m, 
and 500 u.m, respectively. 



4.4. Reconstruction of the SPIRE counts 

We can reconstruct the SPIRE counts using the information 
provided by the S24 + z catalog, the mean color, and the 
standard deviation provided by the stacking analysis. In this 
analysis, we assume that the distribution of the SPIRE flux 
density for a given 24 (im flux density is log-normal (see 
Fig. 4 and Sect. 4.2). For a small scatter (<<1), the standard 
deviation of the logarithm of the flux density cn og _ norm co i or can 
be computed from the standard deviation of the flux density 
o-coior: cr log _ normcolol . = cr color /ln(10). For larger scatter, this 
approximation is no longer valid. However, there is a bijective 
link between the following two pairs of parameters: the mean 
and the scatter of the color in linear units and the same thing in 
logarithmic units. We can thus deduce the two parameters of the 
log-normal distribution (mean and scatter) of the color from the 
linear mean and standard deviation measured by stacking. 

We generate 1000 realizations of the SPIRE counts using the 
following recipe: 

- We take randomly a value of the scatter (see Sect. 4.2). At 
each realization, we used a single value of the scatter for all 
the flux density and redshift bins. 

- In each flux density and redshift bin, we take randomly one 
value of the S spire IS 24 color following the uncertainties (see 
Sect. 4.1). We obtain a relationship between S 24 and the 
color in each redshift slice interpolating between the centers 
of the 24 fj.m bins. 

- We then compute the mean color of each source using the 
previous relationship. 

- For each source, we draw randomly a SPIRE flux density 
from its 24 //m flux, its color and the scatter on it. We assume 
a log-normal distribution. 

- We then compute the counts from the obtained SPIRE flux 
densities. 

The final counts are computed taking the mean and the standard 
deviation of the different realizations. 

Due to the flux density cut of the 24 u.m catalogs, the SPIRE 
simulated catalogs are not complete at the faint end. If there was 
a single color for all objects, the cut of the SPIRE catalog would 
be the SPIRE/24 (a.m color multiplied by the flux density cut at 
24 u.m. Above this limit, the catalog would be complete (sta- 
tistically speaking), and there would be no sources below this 
limit. However, due to the scatter of the colors, this transition 
is smoother. We call the ratio between the reconstructed counts 
(taking into account the 24 u.m selection) and the input counts 
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Figure 7. Completeness of the SPIRE counts reconstructed by stacking 
in a simple case for various values of scatter on the colors. We assume 
power-law counts (dN/dS oc S" 1 ' 5 ), a mean SPIRE/24 u.m ratio, a log- 
normal scatter with different values and a flux density cut at 24 u.m of 
30 pjy. This figure is discussed in Sect. 4.4. 



the completeness. Berta et al. (2011) used the Le Borgne et al. 
(2009) model to estimate the completeness as a function of the 
far-infrared flux density. We have chosen to use a similar, but 
more empirical, method to estimate the completeness and cor- 
rect for it. 

- We generate a mock 24 iim catalog following power-law 
counts with a typical slope in dN/dS oc S~ L5 (Bethermin 
etal. 2010a). 

- We associate a SPIRE flux density with each source of the 
mock catalog using the real colors and scatters measured by 
stacking. The color of each source depends on its 24 (a.m flux 
density and redshift. 

- In each SPIRE flux density bin, we compute the ratio be- 
tween the total number of sources and the number of sources 
which are brighter than the 24 \±m flux density cut. 



The same type of analysis was performed to compute the red- 
shift distribution of S 250,350, or 500 > 6 mJy sources in COSMOS. 
In this case, there is only one flux density bin (> 6 mJy). 



5. Estimation of the statistical uncertainties 

In Sect. 3 and 4, we explained how we derived number counts 
and redshift distributions above and below the confusion limit. 
We also discussed the uncertainties in the corrections applied 
to our measurements. In this section, we explain how the field- 
to-field variance on our measurements is estimated and how we 
combine these uncertainties with the errors on the corrections. 



5. 1 . Sample variance 

Our study is based on only one or two fields depending on the 
flux density regime. The field to field variance cannot thus be 
easily estimated. We have used the same method based on the 
clustering of the sources as in Bethermin et al. (2010a), which is 
briefly described here. 

5.1.1. Principle 

Spatially, sub-mm sources are not Poisson distributed (Blain 
et al. 2004; Farrah et al. 2006; Cooray et al. 2010; Magliocchetti 
et al. 2011), but clustered. The uncertainty, cr N , on the number 
of sources in a given bin, N, is thus not V/V. In large fields, this 
effect is not negligible, and the clustering of the sources must be 
taken into account (Bethermin et al. 2010a). The uncertainties in 
the clustered case are (Wall & Jenkins 2003) 



cr N 



= ^N* 



+ N, 



with 



y = 



n 2 



(4) 



(5) 



where w(8) is the auto-correlation function (ACF) and £1 the 
solid angle of the field. The effect of the clustering on the 
uncertainties depend only on the field (size and shape) and the 
ACF. 



Several realizations of the colors and scatters are used to 
estimate the uncertainties in this correction. Fig. 7 illustrates 
how the completeness values vary with the scatter of the colors 
in a simplified case, where we assume a single color for all 
the sources (S spire IS 24 = 50) and a 24 |om flux density cut 
of 30 /jjy. As expected, without scatter, the transition happens 
around 1 .5 mJy (50 x 0.03), and the width of the transition 
increases with the scatter. Table B.2 and B.3 provides the 
completeness corrections used in GOODS-N and COSMOS. 
We have cut our analysis in COSMOS at 6 mJy in the three 
SPIRE bands, because the completeness in the higher redshift 
bins is only ~50%. Below this limit, we used the GOODS-N 
field where the 24 /mi catalog is deeper. Following the same 
criterion as in COSMOS, we cut our analysis at 2 mJy in 
the three bands. These cuts are slightly arbitrary, because the 
mean S spire/5 24 color of the sources, and consequently the 
completeness, vary with redshift. Nevertheless, we use the same 
cuts for all redshifts in order to simplify the interpretation and 
the discussion. 



5.1.2. Estimation of the auto-correlation function 

The purpose of this paper is not to study the clustering of sub- 
mm galaxies, but is just to compute, with a reasonable accuracy, 
its effect on uncertainties in the number counts. We measured 
the ACF of the resolved sources for the selection in redshift. 
This measurement is performed with the Landy & Szalay (1993) 
estimator: 



w(9) 



DD-2xDR + RR 
RR ' 



(6) 



where DD is the number of pairs separated by an angle between 
6 — d9 1 2 and 9 + d6/2 in the real catalog, RR the number of pairs 
in a Poisson distributed catalog generated with the mask used 
for the source extraction, and DR the number of pairs coming 
from a source in the real catalog and a source in the random 
catalog. The method used to quickly compute the number of 
pairs is described in Appendix A. 
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We fit our results with the following simple form 
(Magliocchetti et al. 201 1): 



w(0) = A 



1 deg 



i-y 



C 



(7) 



where y is fixed at the standard value of 1.8. This simple form 
does not work at small scales (<2'), where the contribution of 
the clustering between the sources in the same dark matter halo 
is not negligible (e.g. Cooray et al. 2010). We use only the scales 
larger than 2'in our analysis. The integral constraint C is a factor 
taking into account the fact that Landy & Szalay (1993) estima- 
tor is biased for finite size survey. C depends on the size and the 
shape of the field and the value of y, and can be computed from 



i-y 



C = 



JLJLfry **«* 



Q 2 



Combining Eq. 5, 7, and 8, we obtain: 
y = AxC. 



(8) 



(9) 



For our masks, C = 1.72 in the COSMOS field and 7.16 in 
GOODS-N. In order to compute the effect of the clustering 
on our error bars on the number counts, we thus have to es- 
timate A in the various redshift and flux bins used in our analysis. 



5.1.3. Uncertainties in the resolved number counts 

We measure the clustering of the resolved sources from the 
source list produced in Sect. 3. If we use only the source 
in a given flux density and redshift bin, we do not obtain 
sufficient signal to noise. Therefore, we compute the ACF of all 
S spire > 20mJy sources in a single flux density bin, but four 
redshift bins (0<z<0.5, 0.5<z<l, l<z<2, and z>2), and assume 
that the ACF does not evolve too much with flux density. We 
obtained very good fits in each redshift bins at each wavelengths 
(reduced x 2 < 1.3 in all bins). From these fits, we compute 
the value of the y parameter and the sample variance on our 
measurements. The uncertainties coming from the sample 
variance are then combined with the uncertainties coming from 
the correction applied to the counts (see Sect. 3.4). 

Table B.4 summarizes the relative contribution of the 
clustering term (<x clus = -\JyN 2 ) to the total sample variance 

(odus+poi = JcTdus + °"poi = VyN 2 + Nj. This contribution is 
dominant in the low flux density bins (~85%), and decreases in 



brighter flux density bins, where N is smaller. We also compared 
the sample variance with the uncertainties in the corrections. 
This last correction increases the uncertainties by less than 40%. 
We are thus dominated by sample variance. 



5.1.4. Uncertaintes in the number counts measured by 
stacking 

The clustering of the SPIRE sources below the confusion 
limit (< 20mJy) measured by stacking (see Sect. 4) cannot 
be measured directly. In our analysis, we started from the 
24 fan population as a prior. We thus use the clustering of 
this population to compute the effect of the clustering on the 
uncertainties, assuming it is close to the one of the SPIRE faint 
sources. The ACF was measured in the same redshift bins as for 



the resolved sources. We then compute the sample variance, and 
combine it with the uncertainties coming from the completeness 
corrections, the colors, and the scatter. 

Table B.5 and B.6 summarize the relative contribution of the 
clustering to the uncertainties. As for the resolved sources, the 
clustering term dominates the Poisson term in the sample vari- 
ance. In contrast to resolved counts, the errors coming from the 
completeness correction and the uncertainties in the colors and 
the scatter dominate the sample variance. A possible bias, due to 
the assumption that the 24 fim and sub-confusion limit 250//m 
population have similar clustering properties, has therefore only 
a modest impact to our uncertainty budget. 

5.1.5. Uncertainties in the redshift distributions 

The ACF is difficult to measure in small redshift bins, because 
the number of sources is small and the signal-to-noise ratio 
is then poor. For this reason, we have estimated how the 
clustering evolves when we reduce the size of a redshift bin. 
To quantify this effect, we compute the ACF of the 24 ii.m 
catalog (the signal for resolved SPIRE sources only is too low) 
in 1 - dz/2 < z < 1 + dz/2 bins with dz varying from 0.1 
to 1. We find A oc dz' 09 . To compute the uncertainties in the 
redshift distribution, we thus use the ACF measured previously 
to compute the sample variance on the counts in large redshift 
bins (0 < z < 0.5,0.5 < z < 1,1 < z < 2, and z > 2), and 
apply the scaling relation y oc A oc dz' 09 - We then derived 
the sample variance, and combine it with uncertainties in the 
correction factor for resolved counts and the ones coming 
from completeness corrections, colors, and scatter for counts 
measured by stacking. 



6. Validation on simulation 

In order to check the accuracy of our methods used to measure 
the number counts, we have performed an end-to-end simula- 
tion. The clustering of the sources below the confusion limit is 
not well known, and its effect on stacking has been estimated 
in Sect 4.3 with three methods based on the data. We have thus 
chosen to use a simulation with a Poisson distribution of the 
sources, because it is easier to generate. 



6.1. Description of the simulation 

Our simulation is based on Bethermin et al. (2011) model, 
which is a parametric model based on the Lagache et al. (2004) 
spectral energy distribution (SED) library (two populations: 
normal and starburst galaxies). This model uses a simple 
broken power-law evolutionary behavior of the characteristic 
luminosity and density of the luminosity function (LF). The free 
parameters of the model were determined by fitting observed 
counts (including the Herschel resolved counts published by 
Oliver et al. (2010b)), LFs, and the CIB. The model has not been 
modified since the publication of the associated paper. Note that 
this model includes the contribution of strongly-lensed sources 
to the counts. 

A mock catalog, containing the 24 //m, 250 ^m, 350 fim, 
and 500 /im flux densities and the redshift of the sources, was 
generated following the model. We then build a map of the 
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COSMOS field from this mock catalog, the SPIRE noise map, 
and the SPIRE PSE We then redo all the analysis described 
Sect. 3 and 4 using the S24 + z mock catalog and the simulated 
SPIRE maps. The SEDs of this model were not calibrated 
following the correlation between stellar mass and the star 
formation rate (roughly proportional to the infrared luminosity) 
and are thus not valid below 8 //m rest-frame. We thus cannot 
use this simulation at redshifts larger than 2. 



6.2. Results 

The Fig. 8 shows the results of this simulation. The recovered 
counts (triangles and diamonds) nicely reproduce the shape 
of the counts. The flux density regime probed by stacking is 
well reproduced (reduced^- 2 = 1.4 for 81 degrees of freedom). 
Paradoxically, the resolved counts are not as well reproduced, 
with some points deviating at more than 3-cr. At bright flux 
densities, our recovered counts are systematically lower than our 
results. It could be due to the fact that the extraction technique 
shares the flux of a bright source between several prior positions. 



7. Number counts 

7.1. Results 

From the extraction with priors presented in Sect. 3, we build 
number counts per redshift slice down to 20 mJy at all three 
SPIRE wavelengths. Thanks to the stacking of the 24 urn 
sources in the COSMOS and GOODS-N fields, we reach 6 mJy 
and 2 mJy, respectively. We checked that the counts deduced 
from stacking analysis are in agreement with resolved counts 
above 20 mJy, but they have larger uncertainties than the 
resolved ones. The COSMOS and GOODS-N counts deduced 
by stacking analysis are also in agreement where they overlap, 
but with much smaller uncertainties in COSMOS due to the size 
of this field. We thus use GOODS-N points only at faint flux 
densities, which the COSMOS data does not constrain. Fig. 9 
and 10 show our results. The points obtained by stacking in the 
GOODS-N and COSMOS fields disagree at 2 cr in the 1 < z < 2 
bin at all SPIRE wavelengths (see Fig. 10), which could be due 
to field-to-field variance. 
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Figure 8. Validation of our method of measurement of the counts at 
250^m (top), 350pm (center), and 500 yum (bottom) from a simulation 
based on Bethermin et al. (201 1) model. Solid lines: input counts from 
the simulated catalog for various cuts in redshift. Diamonds: resolved 
number counts measured using the same method as for the real data. 
Triangles: number counts measured by stacking using the same method 
as for the real data. 



The depth and small error bars of our counts enable us 
to detect with high significance the peak of the Euclidian 
normalized counts near 15, 10 and 5 mJy at 250 urn, 350 urn, 
and 500 u.m, respectively. This maximum was seen at 250 u.m 
and 350 u.m by Glenn et al. (2010). With our new results, we 
are able to detect this maximum at 500 /an. We also start to 
see a power-law behavior below this peak, which was seen 
previously only up to 160 u.m (Papovich et al. 2004; Bethermin 
et al. 2010a; Berta et al. 2011). Nevertheless, the significance 
of this detection is hard to estimate because of the correlation 
between the points obtained by stacking. 



7.2. Comparison with the previous measurements 

We have compared our total counts with previous measure- 
ments (cf. Fig. 9). At high flux densities (S>20 mJy), our 
counts agree with the counts measured from resolved sources 
of Bethermin et al. (2010b) in BLAST, Oliver et al. (2010b) 
in SPIRE/HerMES SDP fields, and Clements et al. (2010) in 



the SPIRE/H-ATLAS SDP field. Our measurements are also in 
agreement with the stacking analysis of Bethermin et al. (2010b) 
of the BLAST data. Our new stacking analysis of the SPIRE data 
reduces the uncertainties by about a factor 5 compared for the 
BLAST data. Finally, we agree with the P(D) analysis of Glenn 
et al. (2010), except for the 6 mJy points at 250 /vm and 500 urn 
which disagrees by about 2 cr with our measurements. Due to 
the number of points compared (21), we expect to have about 2 
points with 2 cr difference, so this is not significant. 

The good agreement between the counts produced by the 
stacking and the P(D) analysis confirms the accuracy of these 
two statistical methods. It also suggests that the galaxies seen in 
the mid-IR are a good tracer of the sources responsible for the 
sub-mm counts, and justifies a posteriori our choice to use the 
24 jt/m sources as a prior. The mid-IR faint and far-IR bright pop- 
ulation thus constitute a small contribution to the number counts. 
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Flux 

mjy 


All 


Normalized counts (dN/dS 
jyi.5 sr -i 

< z < 0.5 0.5 < z < 1 


x S 2 - 5 ) 

1 <z<2 


z>2 


Stacking (GOODS-N) 


2.1 
3.0 

4.2 


4989 ± 824 

7261 ± 1199 

10405 ± 1337 


699 ± 220 1828 ± 553 

994 ±250 2714 ±941 

1499 ± 307 3856 ± 1035 


1392 ± 445 
2013 ± 542 
2812 ± 598 


1068 ± 356 
1538 ± 442 

2237 ±512 


Stacking (COSMOS) 



6.0 


22037 ± 3228 


3086 ± 1082 


61 16 ±1575 


8918 ± 2273 


3915 ± 1265 


8.4 


25787 ± 2704 


3759 ± 906 


7463 ±1389 


9969 ± 1851 


4594 ± 1064 


11.9 


29044 ±2168 


4574 ± 728 


8834 ± 1284 


10727 ± 1345 


4907 ± 843 


16.8 


31574 ±2527 


5503 ± 723 


9899 ± 1516 


11185 ± 1579 


4985 ± 1035 


Resolved (COSMOS) 


23.8 


23851 ± 1595 


5558 ± 733 


6694 ± 760 


7516 ± 1076 


4082 ±519 


33.6 


20926 ± 1654 


5349 ± 804 


4875 ± 690 


7379 ± 1151 


3323 ± 536 


47.4 


12653 ± 1345 


3878 ± 762 


4036 ± 752 


3114 ±688 


1623 ± 433 


67.0 


8319 ± 1337 


3813 ±928 


2281 ±681 


1777 ±613 


447 ± 289 


94.6 


5780 ± 1431 


3606 ± 1142 


596 ± 432 


1279 ±681 


298 ± 304 


133.7 


528 ± 532 


- 


528 ± 532 


- 


- 


188.8 


1306 ± 1084 


1306 ± 1084 


- 


- 


- 



Table 2. Number counts at 250 iim. The errors take into account the statistical uncertainties, including the clustering effect, and the uncertainties 
in the completeness corrections. For the points measured by stacking, we also take into account the uncertainties in the colors and the scatter. The 
uncertainties in the SPIRE absolute calibration are neglected here. 



Flux 
mJy 


All 


Normalized counts (dN/dS x S 25 ) 
Jy 15 sr-' 
< z < 0.5 0.5 < z < 1 1 < z < 2 


z>2 


Stacking (GOODS-N) 


2.1 
3.0 

4.2 


4709 ± 1342 
6949 ± 1167 
9964 ± 1396 


510 ± 159 1695 ± 572 

761 ± 234 2598 ± 870 

1115 ±363 3802 ±1016 


1341 ± 450 
1912 ±481 

2587 ± 520 


1162 ± 1117 
1676 ± 564 
2459 ± 715 


Stacking (COSMOS) 


6.0 
8.4 
11.9 

16.8 


21510 ±3858 
23820 ±3174 
24402 ± 2274 
24229 ±3158 


2143 ± 525 5083 ± 1309 
2458 ±431 5715 ±996 
2638 ±462 6175 ± 1001 
2579 ±790 6115 ±1648 


8691 ± 2354 
9627 ± 1997 
9875 ± 1445 
9953 ± 2208 


5592 ±2711 
6020 ±2216 
5713 ± 1366 
5581 ± 1325 


Resolved (COSMOS) 



23.8 18652 ±1605 1967 ±434 3660 ± 643 8442 ±1215 4581 ±703 

33.6 15285 ±1448 1927 ±484 3600 ± 708 4995 ± 838 4760 ±811 

47.4 9092 ± 1187 927 ± 346 1606 ± 467 3929 ± 828 2628 ± 620 

67.0 3487 ±828 728 ± 390 327 ± 234 1527 ± 553 904 ±416 

94.6 1163 ±630 354 ±351 98 ± 160 710 ±497 

133.7 170 ±273 - - - 170 ± 273 

Table 3. Number counts at 350 (im. 



7.3. Comparison with the models 

In Fig. 10, we compare our results with a set of recent (>2009) 
evolutionary models: 



- The Bethermin et al. (2011) model was presented in 
Sect. 6.1. 

- The Marsden et al. (201 1) model is also a parametric model 
similar to the Bethermin et al. (2011) one, but using a dif- 
ferent SED library, and taking into account the scatter in the 
temperature of the cold dust in the different galaxies. 

- Le Borgne et al. (2009) carried out a non-parametric inver- 
sion of the counts assuming a single population (Chary & 



Elbaz 2001) to determine the evolution of the luminosity 
function with redshift. 

- The Valiante et al. (2009) model used a large library of star- 
burst and AGNs templates. This model takes into account the 
scatter in the temperature of the sources. The parameters of 
the model were tuned manually. 

- The Gruppioni et al. (201 1) model uses 5 separately evolving 
populations, including 3 populations of AGN. 

- The Rahmati & van der Werf (2011) model is based on a 
modified Dale & Helou (2002) library. This model takes into 
account the scatter in the temperature of the sources. It was 
fit to the 850 |a.m counts and redshift distribution. 
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Flux 
mjy 


All 


Normalized counts (dN/dS 
jyi.5 sr -i 

< z < 0.5 0.5 < z < 1 


xS 2 - 5 ) 
1 <z<2 


z>2 


Stacking (GOODS-N) 


2.1 

3.0 

4.2 


3465 ± 864 

5216 ±783 

7244 ± 1089 


368 ± 178 1262 ± 394 
493 ± 244 1966 ± 500 
726 ±360 2823 ±631 


836 ±217 
1220 ±351 
1592 ± 529 


998 ±716 
1535 ±424 
2102 ±615 


Stacking (COSMOS) 



6.0 


12170 ± 1764 


750 ± 263 


2381 ± 562 


4444 ±817 


4594 ± 1435 


8.4 


11446 ± 1716 


596 ± 322 


2107 ± 783 


4481 ± 872 


4260 ± 1210 


11.9 


9917 ±2089 


465 ±380 


1586 ± 1000 


3976 ± 1363 


3888 ± 1167 


16.8 


7540 ± 2665 


358 ± 443 


1055 ± 1192 


3003 ± 1817 


3122 ± 1478 


Resolved (COSMOS) 


23.8 


6298 ± 675 


602 ± 220 


1023 ± 278 


2258 ± 343 


2413 ± 460 


33.6 


4548 ± 656 


248 ± 138 


483 ± 191 


1493 ± 329 


2322 ±516 


47.4 


1143 ±343 


- 


130± 117 


549 ± 238 


463 ±218 


67.0 


343 ±251 


- 


- 


182 ± 185 


160 ± 169 


94.6 


202 ± 230 


- 


- 


100 ± 162 


101 ± 163 



Table 4. Number counts at 500 (im. 



Note that the Bethermin et al. (201 1), Gruppioni et al. (201 1) 
and Rahmati & van der Werf (201 1) models were already tuned 
using recent Herschel data, including the GOODS-N observa- 
tions used here. 

None of these models manages to fully reproduce our 
measurements. The Bethermin et al. (2011), Gruppioni et al. 
(2011) and Rahmati & van der Werf (2011) models are close 
to the data, and broadly reproduce the shape of the counts, but 
still deviate from the measurements by 3-cr. The Le Borgne 
et al. (2009) and Valiante et al. (2009) models underestimate the 
contribution of z < 1 sources to the counts. The Marsden et al. 
(2011) model overestimates the counts at high z (z > 1). Not 
surprisingly, models which use the most recent Herschel data 
and the redshift-dependent observables (redshift distributions, 
luminosity functions, etc.) provide the best match to our find- 
ings. 



concentrated in a 0.7° xO. 7° region, corresponding to a physical 
size of about 20 Mpc. It could be linked with the three candidate 
clusters of galaxies at z ~ 1.8 found by Chiaberge et al. (2010) 
in the same field. Nevertheless, this overdensity could also be an 
artifact of the photometric redshifts. An effect of the polycyclic 
aromatic hydrocarbon (PAH) features redshifting into the 24 (j.m 
band is possible althought less likely, because this should affect 
neighboring redshift bins, due to the band width of Spitzer at 
24 iim {XI AA ~ 3). 

We also used the stacking analysis presented in Sect. 4.4 to 
estimate the redshift distribution of the Sspire >6mJy sources 
in COSMOS. We find a smaller relative contribution of z < 1 
sources than for the 20 mjy flux density cut at 250 /im and 
350 //m. The behavior at z > 1 is similar to that found for the 
20 mjy flux density cut. 



8. Redshift distributions 

8.1. Results 

From the brighter sources extracted using the 24 \±m prior 
(see Sect. 3), we have built the redshift distribution of the 
sources brighter than 20 mjy at 250 \xm, 350 \xm, and 500 ii.m 
(see Fig. 1 1 and Table 5). We find that the distribution of the 
resolved 250 \xm sources is almost flat up to z ~ 1 and decreases 
significantly at higher redshift. At 350 [J.m, the distribution 
peaks near z ~ 1, and the distribution is flatter at high redshift. 
At 500 /im, the contribution of z < 1.5 sources is smaller 
than at shorter wavelengths. Between z - 1.5 and z = 3, our 
measurements are compatible with a flat distribution, however 
the uncertainties are very large. 

At 250 \xm and 350 //m, we clearly see an excess in the 
0.2 < z < 0.4 and 1.8 < z < 2.0 bins. The structure at z - 0.3 
in COSMOS is well known (Scoville et al. 2007). The excess 
near z - 1 .9 could also be explained by a large-scale structure. 
Fig. 12 shows the position of the sources in a thin redshift 
slice between z — 1.85 and z — 1.9. The sources are strongly 



8.2. Comparison with other measurements 

Chapin et al. (201 1) studied the redshift distribution of isolated 
BLAST sources. Their redshift distributions cannot be normal- 
ized by the surface area (because of the isolation ctriterion), and 
the flux density cuts are different; nevertheless, the trends of 
their distributions and their evolution from 250 |o.m to 500 [J.m 
agrees with our findings. 

Amblard et al. (2010) also produced a redshift distribution of 
S 350 > 35 mjy sources in H- ATLAS from a Herschel color-color 
diagram and using assumptions about the FIR/sub-mm SED of 
the sources. They found a strong peak at z - 2, in complete 
disagreement with our distribution. This could be due to the 
fact that they required a 3 cr detection at 250 \im and 500 \xm, 
which correspond to a S 250 > 21 mjy and S 500 > 27mJy. The 
3-cr criterion at 500 |a.m tends to select high-redshift sources, 
because of the shape of the SEDs, the flux of low-z sources 
decreases rapidly between 250 //m and 500 /im. The method is 
also strongly dependent on the dust temperatures of the sources 
assumed in their analysis, due to the degeneracy between dust 
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Figure 11. Redshift distribution of Sspire > 20 mJy (upper panel) and 
Sspire > 6 mJy (lower panel) sources at 250 urn (blue), 350 (j.m (green), 
and 500 u.m (red) in the COSMOS field. 



temperature and redshift for thermal sources. 



8.3. Comparison with the models 

We compared the measured redshift distributions with the 
predictions of the same models as in Sect. 7.3 (Fig. 13). 
Again no model manages to reproduce accurately the redshift 
distributions of the bright resolved sources (S > 20 mJy). Note 
however that Gruppioni et al. (2011) model reasonably fits 
the data at 250 fim and 350 fim at z < 2.5. All the models 
without strong lensing predict a strong break in the redshift 
distributions at z~2.5, which is not present in the data. The 
Bethermin et al. (2011) model, which includes strong lensing, 
predicts a more consistent slope, although the normalization 
is not correct. It could be interpreted as a clue that the high 
redshift tail is due to lensed galaxies (see e.g. Vieira et al. 2010; 
Negrello et al. 2010), but the contribution of lensed galaxies in 
the Bethermin et al. (201 1) model is negligible for a flux density 
cut of 20 mJy 7 .The redshift distribution of the faint sources 
(S > 6mJy) are globally better modeled, a broad agreement 
being found with the Bethermin et al. (2011) and Gruppioni 
et al. (201 1) models, which are fitted using the most recent data. 
The strong disagreement with the models at bright flux densities 
suggests that the bright end of the luminosity function and/or 
the SEDs of the brightest objects are not well modeled by the 
current studies. Our measurements therefore provide significant 
new constraints for such models. 



7 Note however that the lensed objects dominates the redshift distri- 
bution at z > 2 for a flux density cut of 100 mJy. 



Figure 12. Spatial distribution of S 2 so >20 mJy sources in COSMOS. 
Black dots: all redshifts. Blue boxes: only sources in the 1.85< z <1,9 
range. 



9. Cosmic infrared background 

9. 1 . Contribution of the 24 [im-selected sources to the CIB 

The differential contribution of the 24 |xm-selected sources to 
the CIB at longer wavelengths as a function of redshift is a 
relatively unbiased measurement and places tight constraints 
on evolution models. Measurements were performed at 70 and 
160 \im in Spitzer data by Jauzac et al. (201 1), and at 250 iim, 
350 |im, and 500 (xm by Viera et al. (in prep.). The latter were 
performed in GOODS-N, on a small area. We performed the 
same analysis in COSMOS, obtaining smaller uncertainties 
and a better resolution in redshift. To compute this overall 
observable, we have estimated the total surface brightness 
due to the S24 > 80//Jy sources in redshift slices. The results 
of this stacking analysis is shown in Fig. 14 and given in Table 6. 

As in Vieira et al. (in prep.), we find a peak near z = 1. 
The relative contribution of the z < 1 sources decreases with 
wavelength, and the contribution of z > 1 sources increases. We 
observed 2 peaks at z ~ 0.3 and z ~ 1.9, probably associated 
with the overdensities discussed in Sect. 8.1. We compared our 
results with the predictions of the three models which can take 
into account the 24 ixm selection among the six previously- 
compared ones. The Bethermin et al. (2011) model broadly 
reproduces our measurements. Nevertheless, it overpredicts by 
3-cr the observed values below z - 1 at 500 |xm . The Valiante 
et al. (2009) model predicts a large bump near 2 = 2, which is 
not seen. A smaller bump is predicted by the Le Borgne et al. 
(2009) model. However, this model tends to underestimate the 
contribution of z < 1 sources and overestimate the contribution 
of z > 1 ones. The large bumps in the CIB contribution predicted 
by the Valiante et al. (2009) and Le Borgne et al. (2009) models 
around z ~ 2, caused by PAH features, are not seen in our data, 
although there is a single elevated point at z ~ 1.8 in all bands. 
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Figure 13. Redshift distribution of the Sspire > 20 mjy (upper panels) and Sspire > 6 mjy (lower panels) sources at 250 [im (left), 350 urn 
(center), and 500 (j.m (right). We overplot the models of Bethermin et al. (201 1) in red, Valiante et al. (2009) in green, Le Borgne et al. (2009) in 
violet, Gruppioni et al. (201 1) in orange, Rahmati & van der Werf (201 1) in light blue and Marsden et al. (201 1) in dark blue. 



We are unsure as to the cause of this observed feature, but note 
that, given the width of the MIPS 24 /mi filter, we would expect 
any significant PAH contribution to affect multiple redshift bins 
instead of a single point. 



9.2. Properties of the CIB 

The contribution of the 24 \im sources to the CIB is an in- 
teresting quantity for models. Nevertheless, we want to have 
constraints on the total contribution of the galaxies to the CIB, 
even if the uncertainties are larger. These constraints can be 
derived by integrating and extrapolating our new SPIRE number 
counts. 



9.2.1 . Estimate of the contribution to the CIB from galaxies 

We integrated our counts for different cuts in flux density 
density, assuming the data points are connected by power-laws. 
The contribution to the CIB of the sources brighter than the 
brightest constrained flux bin is less than 2% (Bethermin et al. 



2010b), and is neglected. We estimated our error bars using 
a Monte Carlo method. We used the distribution of recovered 
values of the CIB to compute the confidence interval. We 
adopted this method down to faintest flux density probed by 
stacking. In order to take into account cosmic variance, we 
combined the statistical uncertainties with the 15% level of the 
large scale fluctuations measured by Planck collaboration XVIII 
(2011). 

We also extrapolated the contribution of the sources fainter 
than the limit of our counts. The typical faint-end slope of the in- 
frared counts 8 lies in a range between -1.45 and -1.65 (Papovich 
et al. 2004; Bethermin et al. 2010a; Berta et al. 2011). This is 
also the case for our input redshift catalog, even if we select 
only a redshift slice. We thus assumed a slope of -1.55 ±0.10 to 
estimate the contribution of the flux density fainter than the limit 
of the stacking analysis. The errors are estimated using a MC 
process, which takes into account the uncertainties in the faint- 
end slope. By integrating our number counts extrapolated down 
to zero flux density, we find a total contribution of the galaxies 



to the CIB of lO.B^nWm^sr 1 , 



6.46+|^ nWnr 



The slope a of the counts is defined by dN/dS ex 5" 
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redshift range 



dN/dz (in 10 4 x gal.sr" 1 ) 
S 250 > 20 mjy S 350 > 20 mjy S 5m > 20 mjy 



0.0 < 
0.2 < 
0.4 < 
0.6 < 
0.8 < 
1.0 < 
1.2 < 
1.4 < 
1.6 < 
1.8 < 
2.0 < 
2.2 < 
2.4 < 
2.6 < 
2.8 < 



z<0.2 
z<0.4 
z<0.6 
z<0.8 
z< 1.0 
z< 1.2 
z< 1.4 
z< 1.6 
z< 1.8 
z<2.0 
z<2.2 
z<2.4 
z<2.6 
z<2.8 
z<3.0 



212±28 

498+57 

323±39 

386±59 

380+58 

358±54 

249+50 

152+31 

99+23 

196+45 

95+24 

120+32 

74+21 

98+29 

51+16 



106+19 

110+18 

96+16 

152+25 
176+32 
251+52 
178+56 
160+56 
101+37 
196+70 

98+37 
116+44 

65+26 
133+51 

81+33 



30+9 
10+4 
32+10 
28+12 
29+12 
41+16 
63+30 
36+16 
17+8 
69+29 
63+26 
45+18 
37+15 
74+27 
30+12 



S250 > 6 mjy S 350 > 6 mJy S 500 > 6 mJy 



0.0 < z < 0.2 
0.2 < z < 0.5 
0.5 < z < 0.8 
0.8 <z< 1.0 
1.0 <z< 1.5 
1.5<z<2.0 
2.0 < z < 3.0 
z>3.0 



667+76 

1713+236 

1865+277 

2346+500 

1493+183 

832+216 

477+88 

11+1 



428+71 

793+176 

1090+213 

1821+417 

1316+262 

807+206 

531+101 

11+1 



145+66 
129+107 
269+184 
707+339 
555+174 
471+149 

314+76 
9+1 



Table 5. Redshift distribution of the SPIRE sources in COSMOS for 
various flux density cuts at the three SPIRE wavelengths. 



redshift range 


d(vl r )/dz (in nW nT 2 


sr 1 ) 




250 urn 


350 yucn 


500 um 


0.0 < z < 0.2 


1.127+0.137 


0.446+0.057 


0.163+0.024 


0.2 < z < 0.4 


2.657+0.295 


1.037+0.123 


0.359+0.048 


0.4 < z < 0.6 


2.279+0.254 


0.935+0.111 


0.305+0.042 


0.6 < z < 0.8 


3.094+0.339 


1.512+0.175 


0.537+0.068 


0.8 <z< 1.0 


3.857+0.421 


2.218+0.253 


0.853+0.104 


1.0 <z< 1.2 


2.612+0.288 


1.602+0.185 


0.622+0.077 


1.2 <z< 1.4 


1.851+0.206 


1.208+0.140 


0.540+0.066 


1.4 <z< 1.6 


1.459+0.165 


1.008+0.118 


0.476+0.060 


1.6 <z< 1.8 


1.088+0.125 


0.748+0.090 


0.376+0.048 


1.8<z<2.0 


1.658+0.187 


1.189+0.139 


0.604+0.074 


2.0 < z < 2.2 


0.871+0.102 


0.664+0.080 


0.375+0.048 


2.2 < z < 2.4 


0.725+0.086 


0.534+0.066 


0.299+0.039 


2.4 < z < 2.6 


0.615+0.075 


0.460+0.058 


0.242+0.033 


2.6 < z < 2.8 


0.689+0.082 


0.531+0.066 


0.282+0.037 


2.8<z<3.0 


0.401+0.051 


0.308+0.041 


0.173+0.024 


3.0 < z < 3.2 


0.207+0.031 


0.162+0.024 


0.093+0.015 


3.2 < z < 3.4 


0.092+0.016 


0.080+0.014 


0.047+0.009 


3.4 < z < 3.6 


0.052+0.011 


0.043+0.009 


0.028+0.006 


3.6<z<3.8 


0.023+0.007 


0.020+0.006 


0.015+0.005 


3.8<z<4.0 


0.018+0.006 


0.015+0.005 


0.007+0.003 



Table 6. Differential contribution of S24 > 80yuJy sources to the CIB as 
a function of redshift. 



and 2.80+^nWrrT 2 sr 1 at 250 \im, 350 |im, and 500 \im, 
respectively. These values agree at better than 1 cr with the 
FIRAS absolute measurements performed by Fixsen et al. 
(1998) and Lagache et al. (2000). 

We estimated the fraction of the CIB resolved into individual 
sources (S > 20 mJy) using our estimation of the total CIB com- 
ing from our extrapolation of the number counts down to zero 



□ Stacking COSMOS 

— Bethermin et al. model 

— Valiante et al. model 

— Le Borgne et al. model 



250 ^m 
S24>80uJy 



I 1 1 1 , T>rirf^t l THfFH 



w? 




Figure 14. Contribution of S24 > 80//Jy sources to the CIB as a func- 
tion of redshift at 250 (J.m (top), 350 (J.m (center), and 500 (j.m (bottom). 
We overplotted the predictions of the Bethermin et al. (2011) (red), 
Valiante et al. (2009) (green), and Le Borgne et al. (2009) (violet) mod- 
els. 



flux density. We found 15%, 11% and 5% at 250 ixm, 350 ixm, 
and 500 \xm, respectively. When we go down to 2 mJy (the limit 
of the stacking analysis), we resolve 73%, 69%, and 55% of the 
CIB, respectively. 

Fig. 15 shows the cumulative contribution to the CIB as a 
function of the flux density cut. We have compared our results 
with the fraction resolved by previous shallower analyses 
(Bethermin et al. 2010b; Oliver et al. 2010b), and find a 1 it 
agreement. 



9.2.2. CIB build-up as a function of redshift 

From our cumulative number counts as a function of redshift 
(see Sect. 7 . 1 ), we can extrapolate the CIB emitted below a given 
redshift, following the methods presented in Sect. 9.2.1. The re- 



18 



Bethermin et al.: Deep number counts at 250//m, 350//m and 500 [im and CIB build-up. 



Resolved (S > 20 mjy) Stacking (5 > 2 mjy) Total extrapolated 



Absolute measurements 



Wavelength 
|i.m 



level fraction 



level fraction 



level Fixsen et al. (1998) Lagache et al. (2000) 



nW m sr 



250 
350 
500 



1 55 +0 - 30 

i. j -0.30 

77+ 016 

"•' -0.16 

14 +0M 

"- 1 -0.03 



15% 

11% 

5% 



7.40 



1.42 
1.43 

4 50 +0<, ° 
1 54 +0M 



73% 
69% 

55% 



10.13!|» 



6.46+ 



C+1.74 

2-80tl 



10.40±2.30 
5.40±1.60 
2.40±0.60 



11.75±2.90 
6.43±1.59 
2.70±0.67 



Table 7. Summary of the contribution to the CIB for various flux density cuts, and comparison with the absolute measurements (which themselves 
have large uncertainties). 
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Figure 15. Cumulative contribution to the CIB as a function of the flux density cut at 250 [im (left), 350 |J.m (center), and 500 \im (right). Red: 
cumulative contribution from our counts. The asterisks represents the fraction resolved at the limit used for our analysis. Cyan: contribution of the 
BLAST sources probed by stacking (Bethermin et al. 2010b). Green: contribution of the sources resolved by SPIRE (Oliver et al. 2010b). Blue: 
contribution of the sources resolved by BLAST (Bethermin et al. 2010b). Violet hatched region: FIRAS absolute measurement of the CIB; a region 
is hatched here if it is in the 1 -cr confidence region of Fixsen et al. (1998) or Lagache et al. (2000). 



suits are presented in Fig. 16 and Table 8. The redshift at which 
half of the CIB is emitted is 1.04, 1.20, and 1.25, at 250 \im, 
350 |xm, and 500 |xm, respectively. For comparison, Le Floc'h 
et al. (2009) measured value z - 1 .08 at 24 ia.m. Berta et al. 
(201 1) performed the same type of measurement, but consider- 
ing only the resolved CIB, and found z = 0.58, z - 0.67, and 
z - 0.73, at 70 \xm, 100 [J.m, and 160 p.m, respectively. As ex- 
pected, the CIB at longer wavelengths is emitted at higher red- 
shift. The predictions of different models are also shown. The 
Marsden et al. (2011) and Valiante et al. (2009) models strongly 
overpredict the contribution of z > 1 sources. The Le Borgne 
et al. (2009) and Rahmati & van der Werf (201 1) models slightly 
underpredict the contribution of z < 2 sources at 350 //m and 
500 fim. The Gruppioni et al. (201 1) model agrees at l<x with the 
measurements, except a 1.5c underprediction at z~l at 500 pm. 
The Bethermin et al. (201 1) models agrees at l<x with this mea- 
surement. Note, however, that it underestimates by lcr the con- 
tribution of z < 2 sources to the CIB at 250 fim and 350 /mi. 



9.3. Spectral energy distribution of the CIB and total 
integrated CIB 

Combining the total extrapolated CIB measured from deep sur- 
veys at various wavelengths, we can produce a fully-empirical 
SED of the CIB (see Fig. 17). We used the values coming 
from resolved counts at 16 |J.m (Teplitz et al. 2011), 24 |xm 
(Bethermin et al. 2010a), 100 \xm and 160 \xm (Berta et al. 
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Table 8. CIB build-up as a function of redshift at 250 [im, 350 |J.m, and 
500 u-m. 



2011), as well as counts measured by stacking analyses at 
70 \xm (Bethermin et al. 2010a), our new results at 250 \xm, 
350 |im, and 500 \xm, and also resolved sources in lensed areas 
at 850 |xm (Zemcov et al. 2010). From these values, we then 
estimate the total CIB integrated between 8 |xm and 1000 |xm: 
27^1 nWm^.sf 1 . We use power-laws to interpolate between 
the data points. To account for the fact that the different data 
points were estimated in similar fields and are thus likely to 
be significantly correlated, we assume a perfect correlation 
between each wavelengths to obtain conservative uncertainties. 



Bethermin et al.: Deep number counts at 250//m, 350/jm and 500 u.m and CIB build-up. 



19 




Gruppioni et al. Bethermin et al. 
Marsden et al. Valiante et al. 
Rahmati et al. Le Borgne et al. 




Figure 16. Cumulative contribution to the CIB as a function of redshift 
at 250 u.m (top), 350 u.m (center), and 500 (im (bottom), and compari- 
son with the models of Bethermin et al. (2011) in red, Valiante et al. 
(2009) in green, Le Borgne et al. (2009) in violet, Gruppioni et al. 
(2011) in orange, Rahmati & van der Werf (2011) in light blue and 
Marsden et al. (201 1) in dark blue. 



Redshift slice 


integrated CIB intensity 




nWirT 2 sr~' 


< z < 0.5 


6.3 ± 1.5 


0.5 < z < 1 


7.9 + 2.2 


1 <z<2 


7.7 ±2.8 


z>2 


4.7 ± 2.0 



Table 9. Contribution of the various redshift slices to the CIB integrated 
between 8/jm and lOOOyum. 



We also derive the contribution to the total CIB from differ- 
ent redshift slices. We use the extrapolated values deduced from 
the counts per redshift slice of Le Floc'h et al. (2009), Berta et al. 
(2011) and our SPIRE measurements. The Berta et al. (2011) 
counts were integrated following the same method as for our 
SPIRE counts. Fig. 17 (colored lines) shows how the CIB SED is 
built up as a function of redshift. The contribution of the various 
redshift slices to the CIB integrated between 8//m and 1000/im 
is given Table 9. 



the following: 6.3 ±1.5 nWm^sr 1 between 
and z - 0.5; 7.9 ± 2.2 nWm^sr 1 between z - u.j wu 
Z = 1; 7.7 ± 2.8 nWm^sf 1 between z — 1 and z — 2; and 



z = 
0.5 and 



z - i; /./+/.» nw m ~sr * 
4.7 ± 2.0 nW irT 2 sr" 1 above z 



= 2. 



10. Discussion 

10.1. Deep source counts in the 250-500 p.m range 

With our new stacking analysis, we have confirmed, with a 
completely independent method, the deep counts produced by 
Glenn et al. (2010) using a P(D) analysis. Unlike for P(D) 
analysis, our stacking approach allows binning in redshift, 
providing new information on the SPIRE sources. 

Our knowledge on the number counts in this wavelength 
interval has dramatically improved in few last years. Before 
BLAST and Herschel, the source counts were very poorly 
constrained by ground-based observations, e.g. 350 |im had 
only three reported ~20mJy sources (Khan et al. 2007). Now, 
thanks to Herschel, they are well constrained between 2mJy 
and 1 Jy. 



10.2. New statistical constraints for the models 

The number counts alone are not sufficient to constrain evolution 
models. In this paper, we have compared different models fit to 
number counts. Some of these models reproduce the number 
counts using incorrect redshift distributions; while here we show 
that in fact, all the models are ruled out by our measurements. 
This highlights how redshift information is crucial in this 
context. The importance of the redshift distributions of the 
sources and of the CIB was also pointed out by Le Floc'h et al. 
(2009), Jauzac et al. (2011), Bethermin et al. (2011) and Berta 
et al. (2011) among others. Le Floc'h et al. (2009) measured 
the counts and the redshift distribution at 24 |j.m with Spitzer 
in the COSMOS field. Berta et al. (2011) produced a large 
collection of observables in the PACS bands (70 \xm, 100 \xm, 
and 160 urn). Here, we provide the same type of observables 
in the SPIRE bands, using a stacking analysis to reach a depth 
similar to Berta et al. (2011), despite having a stronger confu- 
sion. The combination of these three datasets will provide very 
stringent constraints for the next generation of evolution models. 



10.3. Origin of the sub-mm part of the cosmic infrared 
background 

Thanks to the depth and the precision of our new measurements, 
we can now study the sub-mm part of the CIB from an empirical 
point of view. As predicted by most of the models, the mean 
redshift of the CIB increases with wavelength (e.g. Lagache 
et al. 2005). This confirms that the CIB in the sub-mm domain 
is dominated by the high-redshift populations. The extrapolation 
of our counts down to zero flux density provides an estimation of 
the sub-mm CIB in agreement with the absolute measurements. 
Our reconstruction of the properties of the SPIRE sources from 
the mid-infrared and optical data can thus explain how the 
sub-mm CIB was emitted. 

In addition, these redshift distributions will help to interpret 
the CIB fluctuations measured by Herschel (Amblard et al. 201 1, 
also see calibration, flux cut and galactic cirrus discussion in 
Planck collaboration XVIII (2011)) and Planck (Planck collab- 
oration XVIII 2011). In fact, PACS and SPIRE redshift distri- 
butions constrain the emissivities of the infrared galaxies as a 
function of redshift, and will help to break degeneracies between 
these emissivities and the mass of the dark matter halos hosting 
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Figure 17. Spectral energy distribution of the CIB. Black filled stars: our total extrapolated CIB at 250 u.m, 350 vim, and 500 (im. Black filled 
squares: total extrapolated CIB from deep number counts at 16 u.m (Teplitz et al. 201 1), 24 u.m and 70u.m (Bethermin et al. 2010a), 100 \im and 
160|J.m (Berta et al. 2011), and 850 urn (Zemcov et al. 2010). Colored solid lines: contribution of the z<0.5 (purple), z < 1 (dark blue), and 
z < 2 (red) sources to the CIB from the counts measured by Le Floc'h et al. (2009) at 24 |im, Berta et al. (201 1) at 70 u.m, 100 [im, and 160 u.m), 
and in this paper at 250 \im, 350 u.m, and 500 u.m. Colored filled stars: our total extrapolated CIB at 250 (im, 350 \im, 500 u.m for various cuts in 
redshift. The colored stars indicate our new points. The dashed lines correspond to the extrapolation of these contributions below 24 (j.m and above 
500 \im. Cyan solid line: Absolute CIB spectrum measured by COBE/FIRAS (Lagache et al. 2000). Green triangles: absolute CIB measurements 
performed by COB.E/DIRBE at 100 u.m, 140 (im, and 240 \im (updated in Dole et al. 2006). Yellow diamond: absolute measurements of Penin 
et al. (201 lb) at 160 u.m with Spitzer/MYPS. Orange arrows: upper limits derived from opacity of the Universe to TeV photons (Mazin & Raue 
2007). The Berta et al. (201 1), Penin et al. (201 lb), and COBE/FIRAS points have been slightly shifted in wavelength for clarity. 



the star-forming galaxies in the fluctuation models (e.g. Planck 
collaboration XVIII 20 1 1 ; Penin et al. 20 1 1 a). 



11. Conclusion 

Thanks to the sensitivity of SPIRE and the high-quality of the 
ancillary data in the GOODS and COSMOS fields, we have de- 
termined new statistical constraints on the sub-mm galaxies. The 
main results of this work are: 

- We produced deep counts (down to 2mJy), which con- 
firm the previous measurements performed by stacking 
(Bethermin et al. 2010b) and P(D) analysis (Glenn et al. 
2010), and significantly reduce the uncertainties in the mea- 
surements. In addition, we provide number counts per red- 
shift slice at these wavelengths. 

- We measured the redshift distribution of the sources below 
the confusion limit using a stacking analysis. 

- We compared our results with the predictions of the most 
recent evolutionary models, which do not manage to accu- 



rately reproduce our new points. These new constraints will 
thus be very useful for building a new generation of models. 
From our source counts, we also derived new estimates of 
the CIB level at 250 |xm, 350 iim, and 500 nm, in agreement 
and with an accuracy competitive with the FIRAS absolute 
measurements. We also derived constraints on the redshift 
distribution of the CIB. 

Finally, combining our results with other work, we have es- 
timated the CIB integrated between 8 |xm and 1000 \xm, pro- 
duced by galaxies, to be 27+ 



'^nWm^sf 1 
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Appendix A: Measurement of the auto-correlation 
function 

To compute the uncertainties in our counts, we measure the auto- 
correlation function (ACF) using a new method, based on the 
stacking of a density map, containing in each pixel the number 
of sources centered on it. If we stack this map at the position of 
the sources, the expected stacked image M(ff) will be (Bavouzet 
2008; Bethermin et al. 2010b) 

M(0) = (l+w(0))xp s , (A.l) 

where ps is the source density (in sources pixel -1 ) and w(0) 
the ACF. Note that this method is not fully accurate because a 
relative error of 10~ 3 on p s affects w(ff) by an absolute error of 
the same amount. 

We can generalize this method to compute the Landy & 
Szalay (1993) estimator (w(9) = (DD -2XDR+RR) /RR, see 
Sect. 5.1.2). In this case, we use two density maps, one for the 
real sources, and one for a simulated catalog, called hereafter 
"real" and "random" maps, respectively. DD is estimated by 
stacking of the real map at the positions of the real catalog, 
DR by stacking of the random map at the position of the real 
sources, and RR by stacking of the random map at the position 
of the random sources. We then compute the Landy & Szalay 
(1993) estimator from the three stacked maps: DD, DR, and 
RR. This provides an estimate for w(9, <f>), where (6, (p) are polar 
coordinates. To reduce the noise, we compute the mean in 
several annuli. 



This method has a computation time proportional to N SOUKes , 
instead of A^ ources for the naive one. Nevertheless, the com- 
putation time is also proportional N 2 ixels . A small number of 
pixels reduces the range of scales which can be probed. We thus 
use successive rebinning of our density maps to accelerate the 
computation of the ACF over a wide range of scales. 



Appendix B: Additional tables 
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Flux 




Correction factor 




mJy 


< z < 0.5 


0.5 < z < 1 


1 <z<2 


z>2 


250/jm 


23.8 


0.89+0.08 


0.79+0.08 


0.81+0.08 


0.86+0.09 


33.6 


0.91+0.10 


0.81+0.10 


0.89+0.09 


0.87+0.13 


47.4 


0.92+0.15 


0.94+0.15 


0.82+0.14 


0.85+0.20 


67.0 


0.96+0.20 


0.96+0.27 


0.92+0.29 


0.94+0.55 


94.6 


1.02+0.30 


0.96+0.65 


1.16+0.60 


0.89+1.04 


133.7 


- 


0.97+1.19 


- 


- 


188.8 


0.96+0.83 


- 


- 


- 


350/jm 


23.8 


0.70+0.10 


0.67+0.09 


0.80+0.07 


0.81+0.07 


33.6 


0.82+0.15 


0.82+0.12 


0.73+0.10 


0.88+0.11 


47.4 


0.79+0.24 


0.79+0.19 


0.93+0.17 


0.85+0.17 


67.0 


0.95+0.44 


0.67+0.37 


0.85+0.26 


0.86+0.35 


94.6 


1.07+1.85 


0.61+2.00 


1.10+0.74 


- 


133.7 


- 


- 


- 


0.62+1.72 


500 yum 


23.8 


0.72+0.18 


0.63+0.14 


0.65+0.11 


0.75+0.10 


33.6 


0.74+0.32 


0.66+0.21 


0.77+0.18 


0.90+0.16 


47.4 


- 


0.76+0.50 


0.73+0.28 


0.63+0.23 


67.0 


- 


- 


0.94+1.08 


0.85+0.86 


94.6 


- 


- 


0.62+2.09 


0.64+1.49 



Table B.l. Correction factor applied to the resolved counts in the various flux density and redshift bins. 



Flux 

mJy < z < 0.25 



0.25 < z < 0.5 



Completeness correction factor 
0.5 <z< 0.75 0.75 <z<l 1 < z < 1.5 



1.5 <z<2 2<z< 3 



z>3 











250/jm 










2.1 
3.0 

4.2 


1.13+0.22 
1.05+0.13 
1.02+0.07 


1.27+0.31 
1.11+0.18 
1.04+0.10 


2.63+0.17 
1.81+0.28 
1.38+0.22 


3.20+0.04 
2.22+0.24 
1.65+0.26 


1.75+0.44 
1.35+0.32 
1.15+0.20 


1.55+0.15 
1.24+0.10 
1.09+0.06 


1.71+0.24 
1.36+0.19 
1.17+0.12 


2.55+0.33 
1.75+0.43 
1.34+0.30 



350/jm 



2.1 
3.0 

4.2 



1.26+0.22 
1.10+0.12 
1.03+0.06 



1.17+0.27 
1.07+0.17 
1.02+0.09 



2.11+0.35 
1.54+0.32 
1.25+0.24 



1.83+0.33 
1.35+0.30 
1.14+0.21 



1.61+0.37 
1.24+0.30 
1.08+0.19 



1.62+0.35 
1.26+0.28 
1.10+0.17 



1.72+0.25 
1.33+0.20 
1.14+0.14 



9.15+6.24 
3.74+0.30 
2.06+0.32 



500 yum 



2.1 


1.34+0.54 


1.59+0.31 


1.40+0.31 


1.04+0.09 


1.40+0.37 


1.06+0.12 


1.18+0.27 


1.55+0.19 


3.0 


1.12+0.38 


1.25+0.22 


1.12+0.16 


1.08+0.04 


1.17+0.27 


1.02+0.06 


1.06+0.17 


1.24+0.13 


4.2 


1.04+0.24 


1.09+0.12 


1.03+0.06 


1.02+0.02 


1.06+0.17 


1.00+0.03 


1.02+0.09 


1.09+0.07 



Table B.2. Completeness correction factor applied to the counts by stacking in GOODS-N in the various flux and redshift bins. 
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!m and 500 u.rr 


t and CIB build-up. 




Flux 
mJy 


< z < 0.25 


0.25 < z < 0.5 


Completeness correction factor 
0.5 <z< 0.75 0.75 <z<l 1 < z < 1.5 


1.5<z<2 


2<z<3 


z>3 


250/jm 


6.0 

8.4 

11.9 

16.8 


1.41±0.29 
1.19+0.17 

1.09±0.09 
1.04±0.04 


1.20+0.25 
1.06+0.14 
1.02+0.07 
1.00+0.03 


1.44+0.30 
1.18+0.22 
1.06+0.13 
1.02+0.07 


1.49+0.28 
1.24+0.18 
1.10+0.11 
1.04+0.05 


1.67+0.21 
1.28+0.22 
1.12+0.15 
1.05+0.09 


1.24+0.36 
1.07+0.26 
1.02+0.17 
1.00+0.10 


1.36+0.20 
1.14+0.14 
1.05+0.08 
1.02+0.04 


2.41+0.07 
1.70+0.27 
1.34+0.24 
1.15+0.17 


350/jm 


6.0 

8.4 

11.9 

16.8 


1.40±0.37 
1.18+0.28 
1.07+0.18 
1.02+0.10 


1.12+0.11 
1.04+0.05 
1.02+0.02 
1.00+0.01 


1.17+0.15 
1.06+0.09 
1.02+0.04 
1.01+0.02 


1.33+0.20 
1.14+0.12 
1.05+0.06 
1.02+0.03 


1.52+0.27 
1.21+0.26 
1.08+0.20 
1.03+0.14 


1.35+0.19 
1.13+0.14 
1.05+0.08 
1.02+0.04 


1.67+0.28 
1.30+0.25 
1.13+0.16 
1.05+0.09 


3.03+0.59 
2.06+0.13 
1.59+0.35 
1.33+0.36 


500 yum 


6.0 

8.4 

11.9 

16.8 


1.11+0.30 
1.04+0.22 
1.02+0.15 
1.00+0.11 


1.05+0.18 
1.01+0.10 
1.00+0.05 
1.00+0.02 


1.23+0.29 
1.12+0.20 
1.06+0.14 
1.03+0.09 


1.20+0.35 
1.08+0.25 
1.03+0.16 
1.01+0.10 


1.25+0.25 
1.12+0.19 
1.06+0.14 
1.03+0.09 


1.14+0.25 
1.05+0.15 
1.01+0.09 
1.00+0.05 


1.30+0.26 
1.11+0.22 
1.04+0.16 
1.01+0.10 


2.22+0.52 
1.51+0.53 
1.20+0.36 
1.08+0.22 



Table B.3. Completeness correction factor applied to the counts by stacking in COSMOS in the various flux and redshift bins. 



Flux 




C"cluS / <3"poi+clus 








^clus+poi/^tot 






mJy 


< z < 0.5 


0.5 < z < 1 1 


<z<2 


z>2 


< z < 0.5 


0.5 < z < 1 1 


<z<2 


z>2 


250/jm 


23.8 


0.86 


0.85 


0.92 


0.77 


0.90 


0.90 


0.95 


0.86 


33.6 


0.79 


0.72 


0.88 


0.64 


0.86 


0.85 


0.91 


0.81 


47.4 


0.65 


0.59 


0.68 


0.41 


0.81 


0.79 


0.83 


0.77 


67.0 


0.54 


0.39 


0.48 


0.18 


0.77 


0.74 


0.76 


0.72 


94.6 


0.44 


0.17 


0.33 


0.11 


0.74 


0.72 


0.70 


0.72 


133.7 


- 


0.12 


- 


- 


- 


0.71 


- 


- 


188.8 


0.17 


- 


- 


- 


0.72 


- 


- 


- 


350/jm 


23.8 


0.87 


0.89 


0.93 


0.88 


0.92 


0.94 


0.95 


0.92 


33.6 


0.80 


0.84 


0.84 


0.83 


0.88 


0.90 


0.91 


0.89 


47.4 


0.58 


0.62 


0.73 


0.65 


0.81 


0.82 


0.83 


0.82 


67.0 


0.44 


0.26 


0.46 


0.36 


0.75 


0.78 


0.77 


0.76 


94.6 


0.25 


0.12 


0.26 


- 


0.70 


0.79 


0.70 


- 


133.7 


- 


- 


- 


0.10 


- 


- 


- 


0.78 


500 yum 


23.8 


0.83 


0.83 


0.70 


0.85 


0.90 


0.91 


0.86 


0.91 


33.6 


0.60 


0.61 


0.52 


0.77 


0.82 


0.84 


0.80 


0.86 


47.4 


- 


0.30 


0.27 


0.39 


- 


0.77 


0.77 


0.81 


67.0 


- 


- 


0.13 


0.19 


- 


- 


0.72 


0.74 


94.6 


- 


- 


0.07 


0.12 


- 


- 


0.78 


0.78 



Table B.4. Columns 2 to 5: relative contribution of the clustering term to the total sample variance (Poisson+clustering) of resolved counts 
(Ccius/Cpoi+cius)- Columns 6 to 9: relative contribution the sample variance term to the total uncertainties in resolved counts (<x cIus+poi /o" tot ). cr tot 
contains both the uncertainties in the corrections and the sample variance. 
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Flux 

mJy 


0< 


z<0.5 


Cell 

0.5 < 


is / ^poi+clus 

z< 1 1 


<z<2 


z>2 


0< 


z<0.5 


^clus+poi / ^"tOt 

0.5 < z < 1 1 


<z<2 


z>2 


250/jm 


2.1 
3.0 

4.2 




0.79 
0.77 
0.75 




0.95 
0.94 
0.93 


0.84 
0.82 
0.80 


0.77 
0.74 
0.72 




0.56 
0.66 
0.74 


0.63 
0.58 
0.68 


0.47 
0.54 
0.64 


0.44 
0.50 
0.60 


350/jm 


2.1 
3.0 

4.2 




0.75 
0.73 
0.70 




0.94 
0.94 
0.93 


0.84 
0.81 
0.78 


0.78 
0.76 
0.73 




0.63 
0.66 
0.65 


0.60 
0.61 
0.71 


0.46 
0.60 
0.73 


0.16 

0.46 
0.53 


500 yum 


2.1 
3.0 

4.2 




0.69 
0.65 
0.62 




0.93 

0.92 
0.91 


0.77 
0.75 
0.70 


0.76 
0.74 
0.71 




0.44 
0.45 
0.46 


0.62 
0.70 
0.74 


0.55 
0.53 
0.50 


0.22 
0.51 
0.51 



Table B.5. Columns 2 to 5: relative contribution of the clustering term to the total sample variance (Poisson+clustering) on counts measured by 
stacking in GOODS-N (<x cllis /crp ( , i+dlls ). Columns 6 to 9: relative contribution the sample variance term to the total uncertainties in resolved counts 
(Cdus+poi/<7"iot)- crtot contains both the uncertainties in the completeness corrections, the uncertainties in the mean colors and the scatter, and the 
sample variance. 

FlUX ^clus/^poi+clus ^clus+poi/^tot 

mJy < z < 0.5 0.5 < z < 1 1 < z < 2 z > 2 < z < 0.5 0.5 < z < 1 1 < z < 2 z>2 

250/jm 



6.0 


0.91 


0.97 


0.95 


0.88 


0.26 


0.46 


0.28 


0.21 


8.4 


0.89 


0.97 


0.93 


0.84 


0.39 


0.64 


0.40 


0.31 


11.9 


0.85 


0.95 


0.90 


0.78 


0.61 


0.83 


0.61 


0.45 


16.8 


0.81 


0.93 


0.85 


0.69 


0.78 


0.80 


0.57 


0.42 


350/jm 


6.0 


0.88 


0.97 


0.95 


0.91 


0.38 


0.46 


0.27 


0.14 


8.4 


0.84 


0.96 


0.93 


0.87 


0.56 


0.69 


0.36 


0.19 


11.9 


0.78 


0.93 


0.89 


0.80 


0.61 


0.75 


0.53 


0.32 


16.8 


0.68 


0.89 


0.84 


0.71 


0.39 


0.47 


0.37 


0.36 


500 yum 


6.0 


0.74 


0.94 


0.91 


0.90 


0.32 


0.52 


0.41 


0.22 


8.4 


0.60 


0.89 


0.87 


0.83 


0.25 


0.34 


0.41 


0.26 


11.9 


0.46 


0.80 


0.79 


0.74 


0.22 


0.23 


0.26 


0.27 


16.8 


0.33 


0.64 


0.65 


0.61 


0.20 


0.16 


0.18 


0.21 



Table B.6. Columns 2 to 5: relative contribution of the clustering term to the total sample variance (Poisson+clustering) on counts measured by 
stacking in COSMOS (cr clus /o"p i +c iiis)- Columns 6 to 9: relative contribution the sample variance term to the total uncertainties in resolved counts 
(Ccius+poi/Viot)- a" tot contains both the uncertainties in the completeness corrections, the uncertainties in the mean colors and the scatter, and the 
sample variance. 



